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Abstract 

Lattice gas automata with collision rules that violate the conditions of semi- 
detailed-balance exhibit algebraic decay of equal time spatial correlations be- 
tween fluctuations of conserved densities. This is shown on the basis of a 
systematic microscopic theory. Analytical expressions for the dominant long 
range behavior of correlation functions are derived using kinetic theory. We 
discuss a model of interacting random walkers with x-y anisotropy whose 
pair correlation function decays as 1/r 2 , and an isotropic fluid- type model 
with momentum correlations decaying as 1/r 2 . The pair correlation func- 
tion for an interacting random walker model with interactions satisfying all 
symmetries of the square lattice is shown to have 1/r 4 density correlations. 
Theoretical predictions for the amplitude of the algebraic tails are compared 
with the results of computer simulations. 
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I. INTRODUCTION 



Closed, isolated physical systems, whose dynamics is described by a Hamiltonian H(T), 
reach for long times a thermodynamic equilibrium state in which each microstate with total 
energy E has equal weight p(T) ~ 5(H(T) — E): the so-called microcanonical ensemble. 
When brought into contact with a heat reservoir, so that the energy is not fixed but fluctuates 
around an average value, the system is described by the canonical distribution p(T) ~ 
e -PH{V) ^ w jth p the inverse temperature. An essential observation is that in both cases the 
equilibrium distribution is completely known in terms of the Hamiltonian, without the need 
to explicitly solve the dynamics generated by H(T). 

The situation is quite different in the case of driven systems, where the dynamics does 
not satisfy the detailed balance condition, and prevents the system from reaching thermal 
equilibrium, e.g. due to an external driving field or due to heat reservoirs at different tem- 
peratures. An example of the latter is a fluid layer heated from above and cooled from 
below, so that a temperature gradient across the layer is maintained. After long times this 
system reaches a non-equilibrium steady state. The corresponding phase space distribution 
can only be determined by explicitly solving the dynamics, e.g. using kinetic theory ||T|. 

It is helpful to study simple models for driven systems to gain insight in the nature 
of non-equilibrium steady states, and to compare theoretical predictions with the result of 
computer simulations. It is in fact simple to define models with stochastic dynamics that 
violate detailed-balance. 

A class of models that has been studied quite extensively in recent years are driven kinetic 
Ising models with Kawasaki-type spin-flip dynamics || and certain particle hopping models 
H . For a recent review see Ref . Q . Computer simulations have revealed algebraic decay of 
the density-density correlation function, i.e. G(r) ~ A/r n for large r, in the stationary state. 
Although the exponent n can be determined from symmetry considerations alone, using a 
Langevin equation approach @,[|], there is currently no theory available that predicts the 
amplitude A of the tail. 

We propose lattice gas automata (LGA's) as an alternative class of simplified models 
that can be used to study the basic properties of non-equilibrium steady states. But more 
importantly, we present a systematic approximate theory for the large distance behavior of 
the correlation function of conserved quantities. Thus we are able to calculate the amplitude 
of the algebraic tails, starting from the microscopic definition of the model. 

In addition to the type of lattice on which particles move and a required set of local 
conservation laws (particle density, momentum density, etc.), an LGA is defined by a set of 
stochastic transition probabilities that define the stochastic collision rules at each node. In 
the context of LGA's there is a distinction between collision rules that satisfy the condition 
of detailed or semi-detailed-balance || , and rules that violate this condition. Semi-detailed- 
balance models reach for long times a completely factorized equilibrium state that is inde- 
pendent of the transition probabilities. However, to study non-equilibrium steady states of 
driven systems one needs to consider models with collision rules that violate semi-detailed 
balance. Such collision rules are incompatible with a factorized state. Strong violation of 
semi-detailed-balance may even lead to spatial instability and pattern formation [[^H]. 

An advantage of LGA's over Ising-type models is that can be used to model nonequilib- 



rium states of fluids as well. In Ref. [ 10 1 it is explained how non-detailed-balance LGA-fluids 
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are to be considered as generalizations of driven diffusive systems. 

Here we exclusively deal with LGA's having only stable modes, so that after long times 
a spatially homogeneous, but correlated equilibrium state is reached. (Note that we use 
the term 'equilibrium state' as a synonym for 'steady state', to emphasize that we consider 
LGA's driven only through strictly local collision rules that are the same for each node and 
applied simultaneously to each node at each time step.) Such LGA's can be interpreted as 
effective models, whose dynamics represents a coarse-grained, mesoscopic description of a 
physical system kept out of thermal equilibrium. Due to their discrete nature LGA's are 
relatively easy to analyze, and studying their behavior will provide insight in the physics of 
non-equilibrium processes. On the other hand, many authors use LGA's lacking detailed- 
balance to model physical phenomena, without analyzing how the lack of detailed-balance 
may affect the validity of their conclusions. It is therefore important to have a fundamental 
understanding of the statistical mechanics of non-detailed-balance LGA's. 

To describe the correlations occurring in the correlated equilibrium state of non-detailed- 
balance LGA's a microscopic description beyond the Boltzmann equation is required. Busse- 
maker, Ernst, and Dufty |TT[ were the first to derive kinetic equations for LGA's at the level 
of pair correlations, by neglecting three point and higher order correlation functions. This 
theory successfully predicts the magnitude of the pair correlations between occupation num- 



bers at the same or at nearby nodes, as was shown in Ref. by numerically evaluating 
the solution to the kinetic equations, and comparing it with simulation results. 

Here we extend the analysis to large distances, and show that all LGA's lacking detailed- 
balance possess spatial correlations between fluctuations of locally conserved quantities that 
decay algebraically for large distances. This is surprising since the collision rules only in- 
volve occupation numbers at the same node: zero range interactions thus lead to infinite 
range correlations. The mechanism that is responsible for the buildup of these long range 
correlations involves the slow evolution of diffusive or hydrodynamic modes at large scales. 
It is the same mechanism that is responsible for the existence of the well-known long time 
tails in hydrodynamic time correlation functions of equilibrium fluids, and the logarithmic 



density dependence of transport coefficients [12 



The organization of the paper is as follows. In section [TJ we recapitulate the kinetic 



equations of Ref. fllT] in terms of excess correlation functions, and obtain an expression for 
the pair correlation function in terms of diffusive or hydrodynamic modes that resembles 
results derived from the phenomenological mode coupling theory. This expression is analyzed 
for interacting random walkers on a square lattice with x-y anisotropy in section |TTT], and 
with the full square lattice symmetry in section |TV|. In section [V| we discuss a fluid-type LGA 
with full triangular lattice symmetry, which exhibits long range momentum correlations. We 
end with a discussion in section [V|. 



II. RING KINETIC THEORY 
A. Basic definitions 

We consider an LGA defined on a d- dimensional lattice of linear size L. The lattice 
has periodic boundary conditions and contains V = L d nodes. In this paper we will only 
use two-dimensional models with d = 2. At each node r there are b channels (r, Cj) for 
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moving particles with velocity Cj {i = 1, • • • , b). We will consider two specific LGA's in this 
paper: (i) a model defined on a square lattice, where Cj are the nearest neighbor vectors 
(cos(z — l)7r/2,sin(i — l)vr/2) with i = l,---,4, and (ii) a model defined on a triangular 
lattice, where Cj = (cos(i — 1)tt/3, sin(i — 1)tt/3) for i = 1, • • • , 6; in addition there may be 
a channel i = for a rest particle with Co = 0. The absence, respectively presence, of a 
particle in channel (r, q) is denoted by boolean occupation numbers Sj(r) = {0, 1}. 

The state of a node r is denoted by s(r) = {sj(r)}. During the collision step of the LGA 
the precollision state s(r) is replaced by a postcollision state cr(r) at all nodes simultaneously, 
according to a stochastic process with transition probabilities A sa > 0. The 2 b x 2 b matrix 
A sa is normalized: 

E4t = i. (i) 

The collision step is followed by a propagation step during which a particle with post- 
collisional velocity Cj is moved from node r to a neighboring node r + Cj. The combined 
collision and propagation steps constitute a time evolution of the entire LGA from time t to 
time t + 1. 

In most LGA's the collision rules satisfy certain local conservation laws. For instance, 
in an LGA describing the diffusive behavior of (interacting) random walkers, the number 
of particles at a node does not change during collision, but the distribution among velocity 
directions does. This conservation law is conveniently formulated in terms of the collisional 
invariant a« = 1. In a fluid-type LGA the local momentum at each node is also conserved 
during collision, and we have di = {1, Cj}. Nonzero transition probabilities A sa > are only 
allowed if 

^0^ = ^20^, (2) 

i i 

or, stated compactly, the matrix A sa must satisfy 

J2 a i( a i ~ s i) A sa = 0. (3) 

i 

The transition matrix A sa is said to satisfy the semi-detailed-balance or Stueckelberg 
condition [O if 

5>- = l- (4) 

s 

The stronger detailed-balance condition, A sa = A as , implies semi-detailed-balance on account 
of the normalization (Q). It can be shown that if Eq. (|j) holds, the equilibrium distribution 
is completely factorized over all bV channels (r, Cj), and only depends on the microscopic 
state through global invariants, like the total number of particles or the total momentum 
0. Since the collision step does not change the value of these invariants, it follows that the 
equilibrium distribution is invariant under the collision step. 
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B. Simple versus repeated ring approximation 



We restrict ourselves in this paper to properties of spatially homogeneous equilibrium 
states in LGA's lacking detailed-balance. The quantities of interest are the average occu- 
pation number or single particle distribution function = (sj(r)) and the pair correlation 
function Qij(r) = (5si(r)Ssj(0)) with <5sj(r) = Sj(r) — fa. 



We give a short summary of necessary results derived in Ref. HTT| . In semi- detailed- 
balance models with zero range interactions the average occupation number equals the Fermi 
distribution, and the pair correlation function has the diagonal form 

^(r) = ^(r)=M(r,0)/ i (l-/ i ), (5) 

showing the absence of spatial or velocity correlations. 

Next we consider models that violate semi-detailed-balance. It is convenient to introduce 
the excess pair correlation function 

C«(r) = ^(r)-Sj(r) J (6) 

where a special role is played by the on-node correlations Cy ee Cy(0). In the so-called 
simple ring approximation the average equilibrium occupations, {/«}, are the solution to the 
stationary nonlinear Boltzmann equation, 

n]'°(f) = - Sl )A S(7 F(s) = 0, (7) 

sa 

where the nonlinear Boltzmann operator Q]'°(f) depends on the average occupations /$ 
through the factorized distribution F(s), defined as 

F(s)=Uft(l-f t ) 1 ^. (8) 

i 

The source of all spatial correlations is the matrix E, which in simple ring approximation is 
given by 

E i;i = tff ee XX^o-j - Ss^A^Fis). (9) 

sa 

Once Eij is known, the on-node correlations Cy can be calculated from the stationary ring 
equation, 

Cij = Rij^iEkt- (10) 



The explicit form of the ring operator, R, given in Ref. |TT], is not needed here. 

At a more sophisticated level, the repeated ring approximation, {f{\ and {Cij} are ob- 
tained as the solution to the stationary (generalized) Boltzmann equation, 

n^(/) + EnS(/)c« = o, (ii) 
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where the term containing = d 2 ^' / dfkdfy describes corrections to fij' . The on- 

node excess correlation function C couples the generalized Boltzmann equation ([TT|) to the 
stationary ring kinetic equation (|I0|), where the source matrix E is now given by 

E*j = ^if + ^ij 2 kiCke + XX 1 ~~ ^ij^eCu, (12) 

with fi| 2 fc£ = d 2 Q 2 f/df k dfe and 1^ = S ik 5 je . Furthermore, = (1 + fi) ifc (U + 

where IL is the unit matrix with (l)y = Sij, and Q is the linearized Boltzmann collision 
operator, defined as, 

"•■f-^-^'A (13) 

For a derivation of these equations as well as a detailed discussion of how to obtain a 
(numerical) solution, we refer to Ref. [11]]. Using the definition in Eq. (§) or (12) it can be 



shown that E satisfies all local conservation laws of the model, i.e. 



(aa\E) = (a\E\a) = did j Eg = 0. (14) 



We have found ||TT] that for models with local conservation laws the numerical difference 
between the simple and repeated ring value for dj is on the order of 10%. Corrections to 
the Boltzmann value for fa, as obtained from Eq. (ffl|) , are even smaller — typically 1%. 



As shown in Ref. [ 1 1| and the next subsection, all spatial correlations in the system are 
linear in the source term E. If this term vanishes, all correlations in the system vanish, and 
the equilibrium state is completely factorized. That the source term does indeed vanish if the 
collision rules satisfy the detailed-balance condition can be seen as follows. As noted above 
Eq. (|5|), fi is a Fermi distribution. Then the single node distribution (^j) satisfies the relation 
F(s)A sa = F(a)A sa . Using normalization ([!]) and semi-detailed balance (|) it follows that 
f^j = 0, and consequently Ey = 0, both in simple and repeated ring approximation. If the 
transition rates do not obey the semi-detailed-balance condition (|j) then f2 2 j° is in general 
nonvanishing. 



C. Mode coupling formula 

Here we are concerned with correlation functions of conserved (hydrodynamic) densities. 
In diffusive models the only conserved density is the number density p(r) = SjSj(r). In 
fluid-type models the momentum density g(r) = J2i c iSi(r) is conserved as well. We denote 
the conserved densities collectively as a(r) = J2i a iSi{r). The hydrodynamic correlation 
functions are then expressed in terms of a scalar product 

G a (r) = (6a(r)6a(p)) = ^^(r) = (aa\g(r)}, (15) 

ij 

where the fluctuation Sa (r) = J2i Qi^M- The Fourier transform of the correlation function 
Oij(r), defined by 
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^■(q) = E^ q - r ^(r) 



(16) 



can be split as 



G ij {q) = 6 ij f i {l-f i )+C ij {q). 



(17) 



The constant contribution on the right hand side comes from the diagonal part defined in 
Eq. (HD, and Cjj(q) denotes the Fourier transform of the excess correlation function defined 
in Eq. In a similar manner, the susceptibility Xa(q) is defined as the Fourier transform 
of QJr), i.e. 



Xaiqj 



-iqr 



and split it into two parts: 



Its diagonal part xt is given by 



X«(0) = Xa + AXa(q). 



xt = EK) 2 / l (i-/,)- 



and the excess part Axa(q) by 

A Xa(q) = J2 a i a Ai( ( l) = (aa|C(q)). 



(18) 

(19) 
(20) 

(21) 



'.i 



The main result of Ref. [|l(]] describes the dominant behavior of the susceptibility at small 
wave number (q — > 0) as 



A x«(q) = E( aQ l^(q)^(-q)) 1 — - x ( q)+^(-q) 0/v(q)Vv(-q)|£), 



(22) 



which has the structure of a mode coupling formula. Here Vv(q) anc ^ ^v(q) are the slow right 
and left (diffusive or hydrodynamic) eigenmodes of the LGA, determined by the eigenvectors 
of the lattice Boltzmann equation: 



e * M (q)+iq.c _ ^ _ Q 
e * M (q)+iq-c _ ]1 _ (] T 



Vv(q) 
^(q) 



0; 
0. 



(23) 



Here fl T is the transpose of the linearized Boltzmann collision operator Q in Eq. fll"3|) . The 
matrices e iqc and 1L are diagonal matrices with elements 5ije t<l ' Ci and respectively. The 
eigenvalue or relaxation rate of the slow mode {Vv^Vvl is ^(q)- For small q it behaves 
as z M (q) ~ q 2 for purely diffusive modes, and as z^(q) ~ q for propagating sound modes. 
The right and left eigenmodes ^(q) and Vv(q) are 6-vectors, or (6+ l)-vectors if the model 
admits states with a rest particle, with components ip^ and ip^. They form a biorthonormal 
set, satisfying the orthogonality relation 



(^(q)|e^ c |^(q)) = ^ ^(q)e^ c ^^(q) = 5, v . (24) 

i 

Note that all inner products in this article are defined without complex conjugation. 
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D. Perturbation Theory 



The small-g behavior of the susceptibility Xa(q) determines the long range behavior of the 
corresponding correlation function Q a (r). We therefore write Ax a (q) as a Taylor expansion 
in powers of the wave number q — |q|: 

A Xa (q) = A X (°)(q) + g 2 A x i 2 )(q) + g 4 A x i 4 )(q) + ■ • • . (25) 

Explicit expressions for the functions Ax^ m )(q) occurring in Eq. (^) can be obtained by 
expanding the eigenvectors and eigenvalues in Eq. (^) in powers of q: 

M<ti = Vf(q) + (iq)^H^) + (^) 2 < 2) (q) + ■ ■ ■ ; 

^(q) = ^° } (q) + WW® + fa)V£ 2) (q) + " " " ; (26) 

*M = 4 0) (q) + (<9)4 1} (q) + te) 2 *S°(q) + ■ ■ • • 

From Eqs. @ and ( 13) it follows that 

E a &ij = 0- (27) 

i 

In other words: the collisional invariants a are left zero eigenvectors of Q. The dimensionality 
of the null space of Q is equal to the number of collisional invariants: one for diffusive models, 
and d + 1 for athermal (without energy conservation) fluid-type models. From Eq. ( |2"T| ) 
we conclude that for q = the left zero eigenvectors of the propagator are Vv(0) = a 
with 2^(0) = 0. These eigenmodes /x, associated with local conservation laws, are called 
slow or hydrodynamic modes. It will be shown below that only pairs \xv of slow modes are 
responsible for singularities (here a discontinuity or anisotropy at q = 0) in the q-dependence 
of the susceptibilities, and hence for the existence of algebraic decay of the pair correlation 
function. 

By expanding Eq. ( p3|) in powers of (iq) we obtain the following hierarchy of equations 
for the left zero eigenvectors: 

= 0; (28a) 

n T ^ = (c e + z^; (28b) 

^ = (q + + [s« + \(c, + ,«)>i°\ (28c) 

where Q T is the transpose of Q, and c« = q • Cj. Similar equations hold for the right zero 
eigenvectors, but with Q T replaced by Q. The bi-orthonormality condition (24) must also 
holds to all powers of (iq), which yields 

¥M 0) ) = V; WfW) + « 0) M^ 0) > + (tfW) = o, (29) 

etcetera. Note that if Qij is symmetric so that Q T = Q, then ^(q) and ^(q) are equal, 
up to a normalization factor. The perturbation equations ( |28a| )-( |2"8c| ) have the general form 
Q T ip^ = 1^ , where the inhomogeneous term I^ 1 ' depends on the unknown eigenvalue zfi' . 
As the matrix Vt T has left zero eigenvectors i/j^, it is required that 

o = (m&w) = WW) (so) 

for all slow modes v. Solving these equations for z^y enables us to determine the eigenvalues 
perturbatively. 
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E. Algebraic Correlations 



Once Ax a (q) is calculated, Fourier inversion of Eq. QT8| ) enables us to calculate the 
spatial correlation functions. In the limit of large system size we can make the continuum 
approximation, 

Ly / dq. (31) 

Here v is the volume of a unit cell in the lattice (v = \\f?> for a triangular lattice; v o = 1 
for a square or cubic lattice), and the q-integration extends over the first Brillouin zone. 
The excess correlation function is then given by 

C ^ ) = WY'L d ^' AxM - (32) 

Combining Eqs. (|25"D and (|32|) we have 

C a (r)=C(°)(r)+C( 2 )(r)+C( 4 )(r) + ---. (33) 
Consider the contribution of the 0(q m ) term in Eq. ( [25]) to C a (r), 

Cf } (r ) = ^_ 1^ rfq ff m e «,r ^ (34) 

If (q) = is isotropic, i.e. continuous at q = 0, then the right hand side of 

Eq. fl3"4|) is essentially a representation of (the m-th derivative of) the Dirac delta-function, 
and therefore all correlations are short-ranged. The situation is very different when Ax^ m )(q) 
is anisotropic, i.e. it depends on q as q — ► 0. A rescaling of q in Eq. (p4[) then shows that 



^ m)(r) " j^y^ L dq rel " Axim)(£l) ' (35) 

where r = r/|r|. Therefore for large r the pair correlation function Q a (r) behaves as 

- ^ ( 36 ) 

with a coefficient A(r) that depends on the direction of r. The value of m is determined by 
the first anisotropic term in the expansion ( p5|) of the susceptibility. The amplitude A(r) 
can be calculated from the microscopic definition of the model by performing the Fourier 
integral in Eq. (|3"5]) . 

In the remainder of this paper we determine m and calculate A{v) for two different 
models: (i) interacting random walkers on the square lattice with an anisotropic transition 
matrix A sa yielding spatial density-density correlations of type 1/r 2 or 1/r 4 , depending on 
whether or not the symmetry between x- and y-directions is broken, and (ii) a fluid-type 
model on a triangular lattice with spatial correlations of type 1/r 2 in the momentum density. 
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III. INTERACTING RANDOM WALKERS WITH X-Y ANISOTROPY 



In this section we discuss an LGA for interacting random walkers on the square lattice. 
The collision rules of the model break the symmetry between the x- and y-direction. We 
choose a model that is still invariant under reflections in both the x- and y-axis, so that no 
average particle drift occurs. Collision rules that break the x-y symmetry are most easily 
formulated in terms of the particle flux J(s) corresponding to a state s, 

J(s) = (37) 

i 

We choose the matrix of transition probabilities as 

A S a = exp[J( S ) • M • J (a)] 8(p(s), p{a)). (38) 

where Z(s) is a normalization constant, 

Z(s) = ^ exp[J( S ) • M • J(<r)] S(p(s),p(a)), (39) 

and M is a diagonal matrix, 

M=(&°V (40) 



If (3 X — (3 y — then the detailed-balance condition is satisfied, and a completely factorized 
equilibrium state exists. For all other choices of f3 x and (3 y — positive or negative — the 
density-density correlations in the correlated equilibrium state decay algebraically. In the 
special case (3 X = (3 y ^ the model has the complete symmetry of the underlying lattice. 
This case will be discussed in the next section. In the remainder of this section we show 
that when (3 X ^ f3 y the correlations are of type 1/r 2 . We derive an analytical expression for 
the amplitude, for the specific interacting random walker model defined by Eqs. (|37|)-([40|). 

The system of interacting random walkers (IRW's) on a (bipartite) square lattice, which 
make a move at each time step, consists in fact of two totally independent subsystems: the 
IRW's initially on the even sublattice C + and those initially on the odd sublattice In 
a single time step all particles on C + move to £_ and vice versa. Therefore equal-time 
correlations can only exist between particles at positions r and r' on the same sublattice. 
Consequently, the difference r — r' always belongs to the even sublattice, so that Q{v) = 
for r G and Q(r) is possibly non- vanishing for r G C + . 

The above features of the bipartite square lattice are contained in the mode coupling 
formula (^) through the existence of two slow modes, both contributing to the long range 
part of the pair correlation function. Let N+(t) and N-(t) denote the total number of 
particles at time t on £ + and £_ respectively, then their difference oscillates in time, and 
Nq = (— Y[N + (t) — iV_(t)] is a conserved quantity, just like the total number of particles 
iV = N + (t) + N-(t). The slow mode corresponding to the conservation of Nq is called 
the staggered diffusive mode; the one corresponding to the conservation of N is the usual 
diffusive mode. 
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The regular diffusive mode has a relaxation rate that for small wave number behaves as 



(see Eq. (|A22|) of App. |A|) 

z D (q) ~ -q 2 z 2 (q) = -(D x q 2 x + D y q 2 y ), (41) 



with diffusion coefficients in the x- and y-direction given by Eq. (|A23|) . To leading order 



the excess susceptibility A%(q) = Ax p (q) contains a contribution from a pair of diffusive 
modes, i.e. 



A x (q) = (ll|^(q)^(-q)) {j—^) (M<dM-<d\&)- (42) 

In App. [A]the left and right diffusive eigenvectors ^-o(q) an d ^(q) are calculated using per- 
turbation theory. For small q the amplitude factors in Eq. (^) are calculated in Eqs. ( |A18| ) 
and (|A27|) with the result 



= 1; ftMq)Vfc(-q) \E) = 2(B x q 2 x + B y q 2 ). (43) 
The factor involving the eigenvalue zd is given by 

1 _ e 2z D i fl ) ~ 2 (D x q 2 x + D y q 2 ). (44) 

In Eq. ( |A25|) of App. [A] the coefficients B a are given explicitly. In the majority of publi- 
cations on driven diffusive systems [|5],|2|,|4| the transport coefficients D a and the coefficients 
B a — which in the phenomenological description represent the noise strength of the fluctu- 
ating force in the Langevin Equation — are simply phenomenological input in the theory. 
In the present paper both sets of coefficients D a and B a are calculated from the microscopic 
definition of the model. 

From Eqs. (f4"2"|)-([4"4D it can be seen that the limit Ax^q) = lim^o A%(q) exists and 
that the dominant part of this contribution to the excess susceptibility is given by 

D x q 2 x + D y q 2 y 

Inverse Fourier transformation yields the contribution Gd{y) of the two diffusive modes to 
the large r behavior of the pair correlation function, 

„ , D x B y -D y B,\ D y g* - D x y* 



However, the staggered slow mode also contributes to the excess susceptibility. It occurs 
at the wave vector 7r = (vr,7r), and is intimately related with the diffusive mode /i = D 
occurring at q = 0. We have r(q+ 7r) = e~ l7T * C T (q) = — T(q), where T(q) = e iq ' c (l + £1) is 
the one-step propagator with eigenvalue e ZD ^\ Then ZD(q + n) = -2r>(q) +m, ^(q+^r) = 
-0£i(q), and -0£>(q + it) = ^(q). It follows that A%(q + 7r) = A%(q). The staggered mode 
at q = 7r gives a contribution to the pair correlation function equal to e~ 47r ' r (?D(r), so that 
the final result for the pair correlation reads 
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G(r) = (1 + e—)Gn(r) = { f D ® % + V Z7 ( 47 ) 

10 x + y odd 

To test the accuracy of our prediction we have performed computer simulations for (3 X = 1 
and fly = 3 at the half-filled lattice, where f% = \ for all four velocity directions. To obtain 
numerical values for the source matrix in the repeated ring approximation, Eq. (0), we 
determined the on-site correlations dj using the methods of Ref . [JTTJ . 

Figure [I] shows a comparison between simulation data and the analytical prediction of 
the large r behavior. The repeated ring theory agrees well with the simulation values over 
the range r G [10, 50]. For large r there is a systematic deviation of the simulation data that 
is a result of the slow diffusive equilibration on large spatial scales, according to r 2 ~ Dt, 
where D is the smallest of the two diffusion constants D x and D v . 



IV. INTERACTING RANDOM WALKERS WITH SQUARE LATTICE 

SYMMETRY 

In this section we discuss the behavior of general diffusive LGA's with collision rules that 
obey all symmetries of the underlying lattice, but violate detailed-balance. An example of 
such an LGA is the model of the previous section in the special case j3 x = j3 y ^ 0. The 
collision rules then obey all symmetries of the square lattice (reflection in x- or y-axis, 
and rotation over multiples of 90°), which implies that second rank tensors are isotropic. 
Therefore the anisotropy giving rise to 1/r 2 correlation for /3 X ^ (5 y , as discussed in the 
previous section, is now absent. However, on the square lattice tensors of rank four contain 
anisotropic parts. In what follows we explain how this anisotropy gives rise to correlations 
decaying as 1/r 4 . 

There are again two slow modes: the usual diffusion mode and the staggered diffusive 
mode. The corresponding eigenvalue zo(q) of the diffusion mode has the form 

z D (q) = -Dq 2 + D 2 (q)q 4 + ---. (48) 

All odd terms vanish because of reflection symmetry. On the square lattice, the so-called 
super-Burnett coefficient D 2 (q) depends on the direction q of the wave vector q, and is 
calculated in Eq. (|A30| ) . It contains an anisotropic term equal to — 2D' 2 q x (jy on account of 



Eq. O). 

To calculate the excess correlation function in Eq. fl42|) we analyze its separate factors. 
The factor containing the eigenvalue 2o(q) behaves for small q as 



1 




(49) 



1 _ e2* D (q) 2Dq 2 { \ D 

The first factor in Eq. ( f42"D equals unity for small q (see Eq. (f43|)). The last one behaves like 

(M<l)M-<l)\E) = 2Bq 2 + 25 2 ( q )g 4 + • • • , (50) 



where the isotropic B and the anisotropic B 2 (q) are calculated in Eqs. ( |A33| ) and ( |A34[ ) ; 
B 2 (ci) contains an anisotropic term —2B' 2 q 2 q 2 . Isotropic terms do not contribute to the 
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algebraic correlations. After collecting terms, the dominant anisotropic contribution of two 
diffusion modes to the susceptibility becomes 




2 



€q 2 y q 2 = -M 2 J 2 y q 2 - (51) 

The large r behavior of the inverse Fourier transform of the anisotropic part of Eq. (]5lD for 
r parallel to the x- or y-axis is given by 

A 



a a 

l + e-^ZIZ- (52) 



The factor (1 + e~ i7r ' r ) accounts for the contribution of the staggered diffusive modes, as 
explained below Eq. fl46|) . This result represents the long range behavior of the pair correla- 
tion function of interacting random walkers with interactions having the full square lattice 
symmetry. The important conclusion of this calculation is that the amplitude of the 1/r 4 - 
tail is nonzero for general choices of (3 X = (3 y ^ 0. Thus the model provides an explicit 
microscopic realization of the scenario that was discussed in the context of the Langevin 
equation by Grinstein et al. ||. The l/r 4 -tail is much weaker than the l/r 2 -tail discussed in 
the previous section, and therefore a comparison with computer simulations would require 
a numerical effort that is beyond the scope of this paper. 



V. FLUID-TYPE MODEL 



In this section we study the spatial correlation functions QapiO) — (ga(?)gp(j)) of the 
momentum densities, with a, f3 — {x, y} in a 7-bit LGA-fluid defined on a triangular lattice, 
that allows for a rest particle state (see subsection IIA) and violates detailed-balance. 

The susceptibility x ttj g(q) is defined as the Fourier transform of the correlation function 
( r ) — (9a( T )9p(0)) with a,P = x,y. We decompose Xa$(<i) into a longitudinal and a 
transverse part, as 

x«/?(q) = q a q/3X£(q) + (<W - qaqs)x±(q), (53) 

where x^(q) = (gi{<£)gi{— q)) and x±(q) — (9±( ( i)9±(~ c l)) are scalar fields with identical 
diagonal parts given by 

Xt = Xl = 3/(1 - /) (54) 

on account of Eq. @. The excess parts A^(q) and Ax_i_(q), given by Eq. ( p2[) with a = q 
and a = c± respectively, are in general different. As we will argue below, the limits for 
q — > of Axeiq) and A^±(q), denoted by Axe and Ax±, are non- vanishing. If Axe ^ Ax± 
then xi°/kq) = li m g ^o X«/?(q) is anisotropic at q = and therefore Q a p{^) ~ 1/r 2 for large 
r. 
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To determine the dominant long range part of G a p(jc) we need to know, according to 
Eq. (1^), the right and left hydrodynamic modes {^(q), ^(q)} and eigenvalues ^(q) of 
the lattice Boltzmann equation, defined through Eq. (p3|) . The collisional invariants in this 
model are a a = {1, c x , c y } or equivalently a a = {a p , a i: a^} = {1, q, c±}, where longitudinal 
(£) and transverse (_L) refer to the q-direction. The set {a a } are the zero left eigenvectors of 
the collision operator Q, and a a are the corresponding zero right eigenvectors, i.e. Q T a a = 
and fla a = 0. The left and right eigenvectors form a bi-orthogonal set, i.e.. (a a \ap) = 5 a p. 

Symmetry properties and the complete set of eigenvectors and eigenvalues are discussed 
in App. [FJ, and summarized in Table [TT], where u = a p , u 2 = and u 3 = a±. It is 
convenient for what follows to show how eigenvalues and eigenvectors transform under the 
inversion q — > — q. We first observe that ^ M (q) = z*{— q) for all modes. Moreover, complex 
conjugation of Eq. (^) shows 

V£(q) = Vw(-q); <(q) = z-«r(q); 

^l(q) = -V»j-(-q); *i(q) = «±(q), (55) 

and the same relations with ip p — > 

For small g the shear mode or transverse momentum mode (// = _L) is 

^i 0) (q) = «±; ^i 0) (q) = a±; zi(q) = -vq 2 , (56) 

and the sound modes (fi = a = ±) are 



^ 0) (q) = a e + cn; s a p ; ^ 0) (q) = K a ± + — a P )', z a (q) = -iqav s - Tq 2 . (57) 



The vectors a Q (a = p, £, _L) are also given in Table ||. The shear viscosity u, the speed of 
sound v s , and the sound damping constant T can be expressed in terms of matrix elements 
Qij of the collision operator, as shown in App. B3, where the higher order coefficients in the 
g-expansion, ipji'\ are al so determined. 

We start with the transverse susceptibilities in Eq. (|22|) and observe that only the pair 
(fiu) = (J-_L) has a nonvanishing overlap for small q, i.e. (c±c±\ip±((i)ip±(— q)) — —1 for 
small q. The excess susceptibility then has the form 

Ax±(q) - -^(^(q)|^|^(-q)). (58) 

Inserting the g-expansion ( |26| ) for ^>j_, and using the relations fl5"5|) and (a Q |£'|a 1 g) = (see 
Eq. flUD), the dominant small-g term in Eq. fl58|) is then 

A X ±(q) = ^ [(rf^ltf } ) - 2^f|E|^ 2) )] . (59) 



These terms are evaluated in Eq. ( |B23a| ) of App. B3 with the result 



Axa(q) = l(^j. (60) 

The eigenvalues — ui^ and of Q and E are calculated in Apps. Bl and B2, and listed in 
Tables and |I| 
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Next consider the longitudinal susceptibility in Eq. fl22"|), where only sound modes (yuz/) 
(a, a') give a nonvanishing contribution for small q, i.e. 



(QQ|Vv(q)?M-q)) 
and the excess susceptibility becomes 



4' 



Ax^(q) = | £ [i?(<7 + a')v s + 2Tq 2 }~ 1 (V» ff (q)|S|^(q)). 



(61) 



(62) 



For parallel sound modes, cr' = a, the denominator yields 2iqav s for small q, and the last 
factor in Eq. (|62"D yields 



(^(q)|E|C(q)) ^ ig f(^ 0) |£|^> + (V^i^S)! = 2iqav s S 10 /u 1 . 



(63) 



The latter equality is derived in Eq. (|B24j ); the coefficient £io, defined in Eq. (P326| ) , is 
calculated in Eq. ( |B33|) in terms of the E^s defined in Eqs. (||) and (|T^). For opposite 
sound modes, a' = —a, the denominator becomes 2rg 2 and the latter factor in Eq. ( |62"D 
yields 



(^(q)|£|C(q)> * q 2 \(^E\^) -2(^\E\^) 



3e4 £ 



+ — 



2£ 



10 



r + ^ 2 (|--; 



1 



(64) 



Here the latter equality is derived in Eq. ( |B30| ) and ( |B31| ), and the coefficient S\\ is calculated 
in Eq. (|B33|) in terms of E^s. Combining Eqs. (0)-(0) yields the final result 



3e 4 E 00 E 00 + 6E W 



8Tuji 



v 2 - 



(65) 



In the case of the simple ring approximation, defined in Eq. @ as Ey = VL^ , the above 
expressions simplify considerably because all diagonal elements are vanishing. It follows 
from Eqs. © and © that 



Ei, 



nff = (1 - 2f t ) - s t )A sa F(s) = (1 - 2fi)Q 



1.0 







(66) 



for all i = 0, 1, • • • , b. We have used the relation (<5crj 



cr; 



fi 



valid for boolean variables <7j. In this case the relevant eigenvalue in Eq. 
€4 = — 2Eis and the excess longitudinal susceptibility becomes 



&Xi 



3-^13 3£^io 



8Ttol ' ATcui 



1 



;i - 2f t )ai + / 2 , 
reduces to 



(67) 



This simplification does not apply in the more general (repeated ring) case. 

In general Xi an d X± are different, unless the collision rules satisfy detailed-balance so 
that E w = En — Ei 3 = 0. By inverse Fourier transformation of Eq. (|53|) we find that the 
asymptotic behavior of the correlation function is given by 
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n / \ (x± - xi)Vs fx 2 -y 



2 ) (68) 



and 



(69) 

An equivalent statement, stressing the isotropy of the correlation functions, is that 

W) = 4^ > ( 7 °) 

G±±{f) = —Qu{f), and Q±e(r) = 0. The labels I and _L here refer to the vector r. In 
Eqs. (§H|)-(0) we may replace xe an d X± D Y ^Xe an d ^x±, respectively, on account of 
Eqs. flU) and (||). 

We performed a computer simulation for the triangular lattice fluid-type LGA defined 
as model III in Fig. 1 of Ref. [|14[], where it was used to study tagged particle diffusion in 
a non-detailed-balance LGA-fluid. Figure |2] shows a comparison between simulation results 
and the theoretical predictions for the amplitude of the algebraic tail. The statistics was 
improved by averaging Qu(r) and G±±(r) over all directions. In particular when the repeated 
ring approximation is used, the agreement is quite satisfactory. 

Although in the limit of long times Qu{r) and Q±±(r) are the same up to an overall sign, 
there is an interesting difference concerning the way in which equilibrium is reached. The 
buildup of Qu{r) is governed by traveling sound modes, for which r ~ t. The buildup of 
G±±{f) however involves the diffusive shear mode, so that r 2 ~ vt. For the data shown in 
Fig. the shear viscosity has the value v ~ 0.2, so that the range over which Q±±(r) has 
equilibrated in T eq = 10 4 time steps is (h>T ec ) 1 l 2 ~ 45 lattice spacings, in agreement with 
Fig. |. 



VI. DISCUSSION 

We have formulated a general ring kinetic theory for lattice gas automata, and used it to 
calculate the pair correlation function for conserved densities. These correlation functions 
have algebraic tails, (?(r) ~ A(r)/r n , for large r. The exponent n can be determined on 
the basis of symmetry considerations alone, using a conceptually simple phenomenological 
Langevin equation approach @||. However, a theoretical estimate for the amplitude A(i) 
can only be obtained by approximately solving the kinetic equations that define the evolution 
in phase space, and analyzing the large r behavior of its stationary solution. This is exactly 
what we did in this paper. 

To test the validity of our approach we performed computer simulations for two different 
two-dimensional models, both violating the condition of semi-detailed-balance. First we 
considered a model of interacting random walkers with different diffusion coefficients in the 
x- and y-direction, exhibiting an algebraic decay of the density-density correlation function, 
G p (r) = (5p(r)5p(0)) ~ 1/r 2 . The second model we considered was a fluid-type model in 
which the correlation function of momemtum densities behaves as Q a p{r:) = (g a (r)gp(0)) ~ 
1/r 2 . In both cases we found good agreement between the simulated and theoretical value 



16 



for the amplitude, in particular when we used the so-called repeated ring approximation, in 
which all pair correlation effects are taken into account in a self-consistent manner. 

Most studies of nonequilibrium states using simple models so far have employed kinetic 
Ising models. Since lattice gas automata are easily implemented and analyzed, as well as 
flexible, they provide an attractive alternative. This holds in particular if one wishes to study 
fluid-type systems in which the momentum density is an additional conserved quantity. The 
algebraic momentum correlations discussed in Section [V] have to our knowledge not been 
observed before, either in computer simulations or in Langevin equation studies. It is an 
interesting question whether such correlations could be detected in real systems, e.g. in 
non-equilibrium states of molecular fluids or of granular media. 

We expect that the techniques used here to analyze lattice gas automata can be extended 
to kinetic Ising models, since the latter constitute just a different class of cellular automata. 
This possibility is under investigation. So far there exists no microscopic theory providing 
the amplitude of algebraic spatial correlations in kinetic Ising models 0. 



APPENDIX A: INTERACTING RANDOM WALKERS 

1. Structure of Q 

The right eigenvectors, u n , the left eigenvectors, u n , and the corresponding eigenvalues, 
— u n of the Boltzmann collision operator Q, are defined by 



VLu n 



-uj„u t 



-0J n U, 



(Al) 



The eigenvectors are constructed solely on the basis of the square lattice symmetry, and are 
given by 



1 

2 2 

c * ~ c v 



;i, 1,1,1); 

1,-1,1,-1' 



c x = (1,0,-1,0); 
Cy = (0,1,0,-1). 



(A2) 



Table | shows how these vectors behave under reflection in the x- and y-axis, respectively. 
There are three invariant subspaces, spanned by {1, (? x — c^}, c x , and c y , respectively. Thus, 
on the basis of square symmetry alone, it can be seen that c x and c y are both left and right 
eigenvectors of Q and E (cf. Ref. ||15|| ). The square symmetry implies that Qij has only six 
independent elements, i.e. 



n 



I fi n fi 12 fii3 fii2 ^ 

^21 ^22 ^21 ^24 

fil3 flia Vlu Vl\2 

\ O2I ^24 ^21 ^22 / 



(A3) 



Number conservation, expressed by (l|fi = 0, or explicitly, 

Vln + 21^21 "H ^13 = ^22 + 21^12 + ^24 



(A4) 
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imposes two more relations, leaving only four independent elements. We easily obtain the 
following biorthonormal set of eigenvectors: 



n i2 cl + OnC?, 



'124 f s '21 

U\ = 1; U\ 



2(n 21 cl - n 12 c 2 y 



1/2 2n 



M2 = n 12 + n 21 ; ^ = 4^-c y „ 

«3 = c^; tt 3 = \c x \ (A5) 

M4 = Cy] U4 = ^'-'J/' 

with eigenvalues given by 

uji = 0; co> 2 = 2(f) 12 + O21); w 3 = O13 — fin; co> 4 = fi 24 - f2 22 . (A6) 

The asymmetry f) 12 7^ fi 21 leads to the mixing between the vectors 1 and c 2 —c 2 . If the model 
is symmetric for interchange between the x- and y-directions, then f) 2 i = ^12 > ^22 = f)ii ; 
and D24 = f)i3- 



This is the relevant set of eigenfunctions for the asymmetric interactions of Section |Tl. 
In Section [TV] the interactions do not break the symmetry between the x- and y-directions, 
and the eigenfunctions and eigenvalues simplify. A table similar to Table | which includes 
the behavior under the symmetry x <-> y, can be constructed for this case. All four vectors 
1, c 2 — c 2 , c x , and c y now span one-dimensional invariant subspaces; so that the eigenvectors 
for the symmetric case are: 

u x = 1; u\= jUi, 



U 2 = c\ - c\\ U 2 = \u 2 ] 

u 3 = c x ; u 3 = |« 3 ; (A7) 

U4 = C y ] U4 = 2^4- 

The corresponding eigenvalues are given by 

loi=0; w 2 = 4f) 12 ; uo 3 = uo 4 = 2(Vl l2 + Vl 13 ) . (A8) 
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2. Structure of E 



The eigenvalues e n of the symmetric source matrix E, defined in Eq. (Q) and (0), are 
defined by 



Ev T . 



(A9) 



In the asymmetric case of section III, the inversion symmetry in the x- and y-axis together 
with the symmetry E^ = Eij imposes the structure 



E 



I E n E u E 13 E l2 ~\ 

E\2 E22 E\2 E24 

E13 E\2 E\\ E\2 

\ E\2 E24 E\2 E22 ) 



(A10) 



The conservation laws (l^ll) = (see Eq. (|14])) imposes one more relation between the 
matrix elements E^. Because of the symmetry E = E T there is no distinction between left 
and right eigenvectors. The symmetry argument of Table [I] of course also holds for E. We 
therefore know that 



v 4 



-y ■ 



are two eigenvectors with eigenvalues 

£3 — En — E13; 64 = E22 — E24- 

The two remaining eigenvalues are the solutions of the quadratic equation 



e 2 — e(Eu + E22 + E13 + E24) + (-^13 + En)(Ei2 + -E24) — 



12 



0. 



(All) 



(A12) 



(A13) 



and the corresponding eigenvectors are linear combinations of c 2 . and c 2 . Their explicit form 
will not be needed in the present paper. 

In the symmetric case of section IV, when there is no difference between x- and y- 
directions, we find that Q and E have the same set of eigenvectors, 



with eigenvalues 



Vl 



ei = 0; 



V-2 



2 _ 2 



V 3 



(2 



-4E 



12; 



e 3 = e 4 = En — E 



13- 



(A14) 



(A15) 



Here ei = follows from the conservation law (1|.E|1) = together with the symmetry 
properties of E. 
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3. Asymmetric interacting random walkers 



Inspection of section III shows that we need to calculate the diffusion coefficients D x and 
D y in Eq. fl4f|) and the projections in Eq. fl43|), which also defines the source terms B x and 
By. To determine these quantities we have to calculate 

V> D (q, c) = + (iq)ipi + {iq?^2- (A16) 
The eigenvalue equations ( p3|) show the symmetry properties 



V>£>(q, c) = Vr>(-q, -c) = ^+(q, c) + V-(q, c); 

z D (ci) = z D (-q)~(iq) 2 D + .--, (A17) 

where ^± has an even (+) or odd (— ) parity in q and c separately. With the help of these 
equations we easily obtain 

B x f x + B y q 2 y = ±(M<l)\E\M-ci)) 

= l(iP 1 \E\<il> 1 )-(ik\E\ip 2 ). (A18) 

Here the relations ipo = 1 and = have been used. 

The solution of Eq. ( |28a| ) is ip = 1 = u±, and similarly ip — U\, where u n and u n are 
defined in App. [A f| . We choose the normalization of ^(q) such that its component parallel 
to u\ is unity to all orders in the perturbation. This implies 

(ux\ip D } = (ui\ipo) = 1 

(«i|Vn) = (n>l). (AI9) 

The solution of the second order equation (28_b), where the inhomogeneous term Ip is a 
linear combination of and u^, then becomes 

0i = ^c £ = - ( ^ + (A20) 

Before we can determine ^ we impose the solubility conditions 

(ui\I ( d ] ) = (fiiC/lVi) + (wiP(q) + |c, 2 ). (A2I) 
We obtain in a straightforward manner 

£>(q) = D x q 2 x + D.gJ = -fel^r + \\c t ). (A22) 
Working this out yields 

= n u / j_ _ a, 
* i) 12 + o 21 U 3 V' 

B » - ?rar - • (A23) 
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The eigenf unction can be solved from Eq. (|28c|) as 

i ( D x q 2 x D y q } 



'H-2 



12 



n 



u 2 . 
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(A24) 



Inserting the results of Eqs. ( |A2(j ) and ( |A24j ) into Eq. (A.18) and using the properties of the 
i?-matrix in App. A2 allows us to obtain the coefficients B a in Eq. ( |A18|) as 



B r = — rrrr—(E u + E n — E 22 — E 2i ); 



cui 2tt 



12 



By = -Aj + — - (En + -&L3 — -^22 — ^24)- 
Co>4 2ii21 

In deriving this result we have used the relations = and 

(l\E\u 2 ) = 2(Eu + E r3 — E 22 — E 2 4). 
The projection in Eq. ( f£f ) simply yields 

(ll|^(q)^ D (-q)) = [(l|^ )] 2 = l, 



(A25) 



(A26) 



(A27) 



because of the normalization ( |A19| ). 



4. Symmetric interacting random walkers 

The required calculations for Section [IV] are much more involved, as the eigenvalue ^(q) 



in Eq. (48) and the source terms in Eq. (|50|) must be evaluated up to terms 0(q i ). This 
requires the perturbation equations for t/jqj ipi, ■ - - , "04 • However, a substantial simplification 
occurs because the linearized Boltzmann collision operator is now symmetric, Q = Q T , and 
right and left eigenfunctions are the same. 

The calculations are tedious but straightforward, and we quote some intermediate as 
well as final results. The coefficients of the diffusion mode ^(q) = X)n=o(*?)'Vn are 

ipo = 1; ipi = — —ct; 

UJ 3 

i>2 = — (4-\); (A28) 

UJ 2 

1 2D 
^3 = (4D0 - ±)4 + — (D + G)c t . 

LU 3 UJ 3 

Here the transport coefficients are 



D =i(£-0< (A29) 



and the super-Burnett coefficient is found as 

D 2 (q) = D' 2 (q A x + q 4 y ) + D' 2 ' = -2D' 2 q 2 x q 2 y + D' 2 " . (A30) 
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It has an anisotropic part D' 2 and an isotropic part, whose explicit form is not needed in 
this paper. Only the former part enters into the coefficient (see Eqs. fl50|) and (|5l|) ) of the 
algebraic tail ~ 1/r 4 of the pair correlation function and it reads 

D£ = 4D(D6-i). (A31) 

The source terms in Eq. (|50"D are determined by 

B = |(Vi|^>; S 2 (q) = |(^2|^|V 2 ) - (V>i|£|V>3>, (A32) 
and the results ( |A28|) enable us to calculate these contributions as 



1 e 3 

B = ^( c i\ E \ c t) = ~S\ 

B 2 (q) = B' 2 (q 4 x + #) + B'> = -2B' 2 q 2 x f y + B' 2 " . (A33) 



Again, only the anisotropic part enters the amplitude of the algebraic tail, and is given by 



B' 2 = ^^ + ^(8DQ-I). (A34) 

LU 2 LU 3 

Finally the projection (|5CiD is again given by Eq. ( [A27Q . 

APPENDIX B: FLUID-TYPE MODEL 

The goal of this appendix is to calculate the left and right hydrodynamic modes (^(o) 
and V'm(q) f° r small wavenumber q) of the lattice Boltzmann equation (p3|) for a fluid- type 
LGA, defined on a triangular lattice, that conserves both particle number and momentum 
during collisions. More specifically, we need to determine the expansion coefficients tfj^ 
(n = 0,1,2) and -ip^ (n = 0) of these modes in powers of iq, as well as the coefficients 

A basic ingredient in this calculation is the structure of the linearized Boltzmann collision 
operator fiy (i = 0,1, ■■•,6) in Eq. ([13]), which is a non- symmetric matrix because the 
LGA under consideration violates the semi-detailed-balance condition (H). The appendix is 



organized as follows. In Section |B 1| the eigenvectors and eigenvalues of f2 are calculated in 
terms of its matrix elements, using the triangular symmetry of the lattice; in Section [B2 
the same is done for the symmetric source matrix defined in Eqs. (§) and ([12]) . Section 



calculates the coefficients ip^ and (ip^ \E\ip^)) in so far as they are needed in the body 
of the paper. 

1. Structure of Q 

The left and right eigenvectors, u n and u n , and the corresponding eigenvalues, — u> n , of 
the non-symmetric Q are defined as 

tt T U n = —U n U n ; VtU n = -UJ n U n , (Bl) 
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where (fl T )ij = Qji and i,j = 0,1,' 
biorthogonal set, normalized as (u n 
lattice impose the general structure 



• • , 6. The left and right eigenvectors together form a 
u m ) = 5 nm . The lattice symmetries of the triangular 



n 



( Cto OL\ OL\ OL\ OL\ OL\ OL\ \ 

a,\ a (3 7 5 7 f3 

a,\ (3 a j3 7 6 7 

&i 7 (3 a (3 7 5 

&i 5 7 (3 a (3 7 

ai 7 <5 7 f3 a (3 

(3 7 5 7 (3 a J 



(B2) 



where the submatrix {f^-; z, j = 1, - ■ ■ , 6} is symmetric. We frequently use a notation where 
a 7- vector t>(c) with components Vi{c) = v (c,) (z = 0, 1, • ■ ■ , b) will be denoted as (v(c )|v(c)) 
or {v (c )|f (cj)). The first component f(c ) refers to the rest particle state with c = 0, and 
the remaining components (i = 1, 2, • • • , 6) refer to moving particle states. The conservation 
law (0) implies that the set of collisional invariants, 



{a p ,a £ ,a ± } = {l,c e ,c±} = {110,112,^}, 



(B3) 



(see Table |T[) are left eigenfunctions, i.e. Q T a a = 0. Multiplication of the matrix ( p2|) on 
the left with u = (1|1) and u 2 = (0|q) imposes the conditions 



«o + 6«i = 0; 
a x + a + 2(3 + 2>y + 5 = 0; 
a + {3 -7 - 5 = 0. 



(B4) 



Because of the symmetry fiy = f^-j for z, j = 1, 2, ■ ■ • , 6 the eigenvectors u n (n = 2, 3) are 
proportional to right zero eigenvectors U2 = &£ and tt 3 = a± (see Table [IT]). To construct the 
right zero eigenvectors Uq we note that by symmetry it must have the structure Uq = (xq\x), 
and satisfy Quq = with (uq\uq) = 1. The result is 



(«i|fii) 



«i — «q 



-(ai| - \a Q ). 



(B5) 



It is possible to identify the components of Uq in terms of the equilibrium distribution 
function /i(p) = (fo(p)\f(p)), which is the stationary solution of the nonlinear Boltzmann 



equation (^), Q { ' (/(p)) = 0, at a given density p = fo + 6/ |IB| . Then we have for an 
infinitesimal change <ip in the density: 



= $>°(f(p + dp)) = £ %(/(p))|^p 



(B6) 



where we have used the definition of f^- given above in Eq. ([13|). This allows us to identify 

(B7) 



u 



<ip dp 
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where the speed of sound v s is defined by 



with p = J2i cjifi = 3f = ~(p — f ). One also verifies from Eq. (|B7|) that (u \u ) = 1. From 
the identification of Eqs. ( |B5[ ) and ( |B7[ ) we obtain 

,2 _ i f Q o \ _ i ( ^oo 

-«o - ai J 2 V^oo - ^01 



= £ h-^- =1 T^V" • (B9) 



Table || shows the complete set of right and left eigenvectors of f2. Most eigenvectors can 
be found on the basis of symmetry arguments alone, except u\ and u\. The corresponding 
eigenvalues are found as 

^1 = ^10-^00; ^4 = ^5 = 2(^12-^14); ^6 = 3(^12-^13)- (BIO) 



2. Structure of E 

So far we have calculated the eigenvectors and eigenvalues of Q. In a similar manner we 
can do so for the source matrix defined in Eqs. (|9]) or (|I2]). Then 

Ev n = e n v n . (Bll) 

As Eij = Eji is symmetric, left and right eigenvectors are the same up to a normalization 
factor, i.e. v n = v n /(v n \v n ). The structure of E is also given by Eq. (p32|) , with a% = a,\. 
The lattice symmetry of the submatrix {Eij]i,j = 1,2, ■••,6} implies that v n = u n for 
n — 2, 3, • - • , 6. The conservation laws (O) imply 



(1|£|1> = ( Ci \E\c e ) = (c ± \E\c ± ) = (B12) 

which imposes two relations between the matrix elements E^, and implies that t>2 = Q and 
^3 = c_i_ are zero-eigenvectors with €2 = £3 = 0. However, eo 7^ 0, and the corresponding 
eigenvector vq is a linear combination of uq and u\. The remaining eigenvalues in terms of 
E^ are obtained from Eqs. ( B1(J[ ) and ( |B12j ) and read 



e 4 = e 5 = 2(_Ei4 — Eu) = 2(En — Ei 3 ) 

e 6 = 3(^x3 - E 12 ) = 3(E n - E u ). (B13) 



The results are summarized in Table 111. 



3. Perturbation theory to 0(q 2 ) 

In this appendix we calculate the hydrodynamic modes and matrix elements occurring 
in Eqs. fl62|) and (p8|) by means of perturbation theory for degenerate eigenvalues. We use 



as a basis the eigenvectors u n and u n of fl that were constructed in subsection B 1 . Our 
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starting point is the observation that according to Eq. ( |28a| ) a hydrodynamic zeroth order 
mode ip^ will be a linear combination of collisional invariants: 

^i 0) = C„ p a p + C^at + C p± a ± = £ C pf3 a p . (B14) 

P 

The solubility condition (|30|) for tM°) yields the equation ^p{joL a \ci + z^\ap)C^ = with 
a — p,£, _L Nonzero solutions C pa only exist if satisfies the secular equation, 

det\(a a \c e + z i p 1) \ap)\ = 0. (B15) 

The only nonvanishing elements are (a p \ci\ai) = v 2 s and (ai\ci\a p ) = 1, and furthermore 
(a a \ap) = 5 a p. From the secular equation we find three non-degenerate eigenvalues and 
subsequently we can determine the coefficients C pa , up to a normalization constant. Thus 
the threefold degeneracy of the null space of Q is lifted, and we have a right shear mode 
fi = ±, 

V>i 0) = c ± z { l ] = 0, (B16) 
and a pair of right sound modes \i = a = ±, 

^ 0) = a £ + <ro s a p = -av s . (B17) 

In a similar manner we obtain the corresponding left eigenvectors of Eq. ( |B15| ) as 

7(0) - i 
if>± = a± = gC^ 

ffi = l(a t + -a p ) (B18) 

Next we consider the first order left hydrodynamic modes. The right modes ipj™^ are only 
required to zeroth order (n — 0). The formal solution of Eq. fl28b[ ) is 

^ = ^(Q + z^f + £ B&tyf . (B19) 

It contains an arbitrary linear combination of zeroth order modes. We always choose the 
normalization such that the projection of ip p onto ipfi) is unity, i.e. (ijj^ \tp p ) = (ip^ |V>1°') = 1 
and consequently for n > 1, 

(?l< n) } = or flg = 0. (B20) 

The corresponding coefficients Bj^J in the right eigenmodes ip^ 1 ' are then determined by 
the normalization conditions (p9|). They are however not needed in the present paper. The 



remaining coefficients Bty (A ^ /i) as well as the next order eigenvalue z^f> are determined 
from the solubility conditions of the second order perturbation equations (PU|). The method 



to determine B^ has been explained in Refs. [T7| for a symmetric collision operator f2, but 
the steps are all very similar. In this manner we find for the eigenvalues to relevant order: 
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-vq 
-iqav s 



rY, 



(B21) 



where v s is the speed of sound, given in Eq. ( |B9| ) , while r = \{y + () is the sound damping 
constant, with v and ( the shear and bulk viscosity, respectively, given in terms of the 
eigenvalues — u n as (see Table [TT]) 



v 



1 



c 



1 



The eigenmodes to first order are found as 





i 

= W 




i 

~ W 




i 







c e c± 



-U4 



OJ4. 

(q - <rv 8 )il>® + B D 

, n 1, err , 



I' 



2v s 



(B22) 



(B23a) 



(B23b) 



«5 Mi r ar 

+ 77 M o - 

Cg> 5 CUi 2 2t> s 



where Table [TT] has been used. The coefficients 5j_a = and B a \ = for A 7^ cr. 

The results so far are sufficient to calculate the transverse susceptibility (|58|)and fl59|). 



From Table |T| we conclude that Eip± = Eu% = 0, and consequently ('^j/'I-^I^'Il'') = 0. The 
remaining term in Eq. fl59p, combined with Eq. ( B23b| ) then yields 



,(°)| 



,(2)\ 



Ax±(q) 



-^(u 4 |w 4 ) 
2z/a; 4 



3e 4 
3 1/0*4 ' 



(B24) 



where _Ew 4 = e 4 w 4 (see Table |T|). 

Next we consider the contributions entering the longitudinal susceptibility in Eq. fl62|), 
starting with the parallel sound modes. To calculate the matrix elements in Eq. ( |63"D we 
observe that Evp^ = av s Euo (see Table plf ). This permits us to combine the terms on the 
right hand side of Eq. (|63| ) into 



(W(q)|£|C CT (q)> = -^(woi^l^ +«> 



05 H 

w 5 



00 



= 2iqav s 

To obtain the second line we have inserted Eq. (P23b|) and introduced 



(B25) 



(B26) 



According to Table [TTJ, n 5 = w 5 * s an eigenvector of .E, orthogonal to the subspace spanned 
by uo and u\\ consequently £05 = 0. Moreover £ 00 = (11\E) = because of the conservation 

To calculate the contributions of the 
Ml). Using E-ip^ = av s Eu we write 



laws (|14]). This yields the result listed in Eq. (^3 
opposite sound modes we need (ip^\E\ip^) in Eq. 



26 



(Vi 0) |J3|V£ 8) > = ctv s (u \E\{uo(u \^) + Wl (%|^ 2) )}) 
= crw s £ i(«i|V'i 2) )- 



(B27) 

m 



To obtain the first equality we note that <?oo = and that all off-diagonal elements £ nn 
(n ^ m), except £ i 7^ (see Table [TTT|) . 

The coefficient (mi|^ 2 ^) in Eq. (|B26| ) can be calculated most conveniently by projecting 
the second order eigenvalue equation ( |28c|) onto the vector u\, i.e. 



(fiil^?*) = («x(q - ™ S M 1} > + r<fi#< 0) ) + |(«x(Q - av s )\{c] - v 2 )). 
Substituting Eq. ( |B23b| ), and using the relations 



(B28) 



u 2 - av s Ui, 



(«i|Vf ) > = 0; 
(u l \n T ^) = -u 1 (u l \^), 



(B29) 



we find for the coefficient 



a 



2v s Ui 



2vl 



(B30) 



The first term in Eq. (|64l) , (tp^\E\ip^) , is obtained similarly by substituting Eq. (B23b|) , 
and yields 



4 u4 W\ Ul 



(B31) 



Substituting Eqs. ( p26|) , ( p30|) , and (p31|) into Eq. (||) then yields the final result, Eq. 
listed in the body of the paper. The excess susceptibility (^) is obtained by substituting 
Eqs. (|63|) and (EM) into (|62|), and carrying out the a-summation with the result 



3e4 



v 2 £ 



10 



16Tu( AFuj( 2rcJi W x 



1 



(B32) 



The final step is to express the coefficients Sn and Sio in Eq. (|B25| ) in terms of the matrix 
elements defined in Eqs. (||) and (|12|). To do so, we write U\ — (| — f 2 )^o — |(1 — c 2 ) 
and use £qo = to find 



£ w = -i(l\E\l - c 2 ) = -|(E 00 + 6£ 10 ); 
£ n = 1(1 - c 2 \E\l - c 2 ) - (| - t; 2 )(l|E|l - c 2 ) 
= iE o-(|-^)(Eoo + 6E 10 ). 



(B33) 
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TABLES 



TABLE I. Symmetries on the square lattice. 
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u n (c) 


X <r+ —x 




1 
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+ 






+ 


2 




c 2 - c 2 
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+ 


3 












+ 


4 




Cy 


+ 










TABLE II. 


Left and right eig 


envectors u n i = u n (ci) with i = 


0,1," 


, b, and eij 


renvalues — w n 


of linearized B< 


altzmann operator f2 for 7-bit fluid-type LGA on triangul. 


ir lattice 




n 


u n (c) 


u n (c) 


-w„ 



1 


a p = 1 

h 2 - 


~a p = (l-2v 2 s \\v 2 s ) 
(— 2|i) 





2 




a e = C£ 


an = \u 2 







3 




a± = c± 


a± = g«3 







4 
5 
6 


C£C± 

r 2 ±r 2 

£ 2 

(o|(-) m ) 


|« 4 
|tt 5 
g«6 


— LL>4 

-w 6 



TABLE III. Eigenvectors v n (c) and eigenvalues e n of source matrix E for 7-bit fluid-type LGA 



on triangular lattice 



n 


V n 


V n 


e n 





linear comb. 


linear comb. 




1 


uq and u\ 


of uq and u\ 


ei 


2 


U2 


U2 





3 


U3 


U 3 





4 


«4 


M4 


€4 


5 


U 5 


U 5 


€5 
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Uq 


Uq 


£6 
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FIGURES 
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FIG. 1. Anisotropic interacting random walker model. Pair correlation function G(r), at even 
r values, with r = (r, 0) along the x-axis, for interacting random walkers on a square lattice with 
interactions that break the symmetry between the x- and y-axis. For r = odd the pair correlation 
function vanishes. The average density per velocity channel is / = ^, and the model parameters 
are f3 x = 1 and (3 y = 3. Symbols with error bars indicate simulation results for a system of 
512 2 nodes, with an equilibration time of T cq = 10 4 time steps. The lines denote the asymptotic 
algebraic tail ~ 1/V 2 , as predicted by ring kinetic theory in the simple (dashed line) and repeated 
ring approximation (solid line). 
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FIG. 2. Isotropic fluid-type model denned in section [V| Correlation function for the longitu- 
dinal, Qu{f)) and transverse, G±±(r), components of the momentum density. The average density 
is / = \ and the total momentum is zero. Symbols with error bars indicated simulation results for 
a system of 256 2 nodes, with an equilibration time of T eq = 10 4 time steps. The lines denote the 
asymptotic algebraic tail ~ 1/V 2 , as predicted by ring kinetic theory in the simple (dashed line) 
and repeated ring approximation (solid line). 
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